﻿ using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using NUnit.Framework;
using System.Diagnostics;

namespace AAPlus.Tests
{
    [TestFixture]
    public class AAPlusTests
    {
        void LunarRaDec(double JD, out double RA, out double Dec)
        {
            double Lambda = Moon.EclipticLongitude(JD);
            double Beta = Moon.EclipticLatitude(JD);
            double Epsilon = Nutation.TrueObliquityOfEcliptic(JD);
            Coordinate2D Lunarcoord = CoordinateTransformation.Ecliptic2Equatorial(Lambda, Beta, Epsilon);

            RA = Lunarcoord.X;
            Dec = Lunarcoord.Y;
        }

        [Test]
        public void NearParabolicObjectElementsTest()
        {
            NearParabolicObjectElements elements6 = new NearParabolicObjectElements();
            elements6.q = 0.921326;
            elements6.e = 1;
            elements6.i = 0; //unknown
            elements6.omega = 0; //unknown
            elements6.w = 0; //unknown
            elements6.T = 0;
            elements6.JDEquinox = 0;
            NearParabolicObjectDetails details3 = NearParabolic.Calculate(138.4783, ref elements6);

            elements6.q = 0.1;
            elements6.e = 0.987;
            details3 = NearParabolic.Calculate(254.9, ref elements6);

            elements6.q = 0.123456;
            elements6.e = 0.99997;
            details3 = NearParabolic.Calculate(-30.47, ref elements6);

            elements6.q = 3.363943;
            elements6.e = 1.05731;
            details3 = NearParabolic.Calculate(1237.1, ref elements6);

            elements6.q = 0.5871018;
            elements6.e = 0.9672746;
            details3 = NearParabolic.Calculate(20, ref elements6);

            details3 = NearParabolic.Calculate(0, ref elements6);

            EclipticalElementDetails ed5 = EclipticalElements.Calculate(131.5856, 242.6797, 138.6637, 2433282.4235,
                                                                        2448188.500000 + 0.6954 - 63.6954);
            EclipticalElementDetails ed6 = EclipticalElements.Calculate(131.5856, 242.6797, 138.6637, 2433282.4235,
                                                                        2433282.4235);
            EclipticalElementDetails ed7 = EclipticalElements.FK4B1950ToFK5J2000(131.5856, 242.6797, 138.6637);

            elements6.q = 0.93858;
            elements6.e = 1.000270;
            elements6.i = ed5.i;
            elements6.omega = ed5.omega;
            elements6.w = ed5.w;
            elements6.T = 2448188.500000 + 0.6954;
            elements6.JDEquinox = elements6.T;
            NearParabolicObjectDetails details4 = NearParabolic.Calculate(elements6.T - 63.6954, ref elements6);
        }

        [Test]
        public void DatesOtherTest()
        {
            CalendarDate JulianDate = MoslemCalendar.MoslemToJulian(1421, 1, 1);
            CalendarDate GregorianDate = Date.JulianToGregorian(JulianDate.Year, JulianDate.Month, JulianDate.Day);
            CalendarDate JulianDate2 = Date.GregorianToJulian(GregorianDate.Year, GregorianDate.Month,
                                                              GregorianDate.Day);
            CalendarDate MoslemDate = MoslemCalendar.JulianToMoslem(JulianDate2.Year, JulianDate2.Month,
                                                                    JulianDate2.Day);
            bool bLeap = MoslemCalendar.IsLeap(1421);

            MoslemDate = MoslemCalendar.JulianToMoslem(2006, 12, 31);
            CalendarDate OriginalMoslemDate = MoslemCalendar.MoslemToJulian(MoslemDate.Year, MoslemDate.Month,
                                                                            MoslemDate.Day);
            MoslemDate = MoslemCalendar.JulianToMoslem(2007, 1, 1);
            OriginalMoslemDate = MoslemCalendar.MoslemToJulian(MoslemDate.Year, MoslemDate.Month, MoslemDate.Day);

            CalendarDate JulianDate3 = Date.GregorianToJulian(1991, 8, 13);
            CalendarDate MoslemDate2 = MoslemCalendar.JulianToMoslem(JulianDate3.Year, JulianDate3.Month,
                                                                     JulianDate3.Day);
            CalendarDate JulianDate4 = MoslemCalendar.MoslemToJulian(MoslemDate2.Year, MoslemDate2.Month,
                                                                     MoslemDate2.Day);
            CalendarDate GregorianDate2 = Date.JulianToGregorian(JulianDate4.Year, JulianDate4.Month,
                                                                 JulianDate4.Day);

            CalendarDate JewishDate = JewishCalendar.DateOfPesach(1990);
            bLeap = JewishCalendar.IsLeap(JewishDate.Year);
            bLeap = JewishCalendar.IsLeap(5751);
            long DaysInJewishYear = JewishCalendar.DaysInYear(JewishDate.Year);
            DaysInJewishYear = JewishCalendar.DaysInYear(5751);
        }

        [Test]
        public void EclipsesTest()
        {
            SolarEclipseDetails EclipseDetails = Eclipses.CalculateSolar(-82);
            SolarEclipseDetails EclipseDetails2 = Eclipses.CalculateSolar(118);
            LunarEclipseDetails EclipseDetails3 = Eclipses.CalculateLunar(-328.5);
            LunarEclipseDetails EclipseDetails4 = Eclipses.CalculateLunar(-30.5); //No lunar eclipse
            EclipseDetails4 = Eclipses.CalculateLunar(-29.5); //No lunar eclipse
            EclipseDetails4 = Eclipses.CalculateLunar(-28.5); //Aha, found you!
        }

        [Test]
        public void MoonPropertiesTest()
        {
            double ApproxK = MoonPhases.K(1977.13);
            double NewMoonJD = MoonPhases.TruePhase(-283);

            double ApproxK2 = MoonPhases.K(2044);
            double LastQuarterJD = MoonPhases.TruePhase(544.75);

            double MoonDeclinationK = MoonMaxDeclinations.K(1988.95);

            double MoonNorthDec = MoonMaxDeclinations.TrueGreatestDeclination(-148, true);
            double MoonNorthDecValue = MoonMaxDeclinations.TrueGreatestDeclinationValue(-148, true);

            double MoonSouthDec = MoonMaxDeclinations.TrueGreatestDeclination(659, false);
            double MoonSouthDecValue = MoonMaxDeclinations.TrueGreatestDeclinationValue(659, false);

            double MoonNorthDec2 = MoonMaxDeclinations.TrueGreatestDeclination(-26788, true);
            double MoonNorthDecValue2 = MoonMaxDeclinations.TrueGreatestDeclinationValue(-26788, true);
        }

        [Test]
        public void PhysicalObjectTest()
        {
            PhysicalSunDetails psd = PhysicalSun.Calculate(2448908.50068);
            double JED = PhysicalSun.TimeOfStartOfRotation(1699);

            PhysicalMarsDetails MarsDetails = PhysicalMars.Calculate(2448935.500683);

            PhysicalJupiterDetails JupiterDetails = PhysicalJupiter.Calculate(2448972.50068);

            //The example as given in the book
            GalileanMoonsDetails GalileanDetails = GalileanMoons.Calculate(2448972.50068);

            //Calculate the Eclipse Disappearance of Satellite 1 on February 1 2004 at 13:32 UCT
            double JD = 2453037.05903;
            int i;
            for (i = 0; i < 10; i++)
            {
                GalileanMoonsDetails GalileanDetails1 = GalileanMoons.Calculate(JD);
                JD += (1.0 / 1440);
            }

            //Calculate the Shadow Egress of Satellite 1 on February 2  2004 at 13:07 UT
            JD = 2453038.04236;
            for (i = 0; i < 10; i++)
            {
                GalileanMoonsDetails GalileanDetails1 = GalileanMoons.Calculate(JD);
                JD += (1.0 / 1440);
            }

            //Calculate the Shadow Ingress of Satellite 4 on February 6 2004 at 22:59 UCT
            JD = 2453042.45486;
            for (i = 0; i < 10; i++)
            {
                GalileanMoonsDetails GalileanDetails1 = GalileanMoons.Calculate(JD);
                JD += (1.0 / 1440);
            }

            //Calculate the Shadow Egress of Satellite 4 on February 7 2004 at 2:41 UCT
            JD = 2453042.61042;
            for (i = 0; i < 10; i++)
            {
                GalileanMoonsDetails GalileanDetails1 = GalileanMoons.Calculate(JD);
                JD += (1.0 / 1440);
            }

            //Calculate the Transit Ingress of Satellite 4 on February 7 2004 at 5:07 UCT
            JD = 2453042.71181;
            for (i = 0; i < 10; i++)
            {
                GalileanMoonsDetails GalileanDetails1 = GalileanMoons.Calculate(JD);
                JD += (1.0 / 1440);
            }

            //Calculate the Transit Egress of Satellite 4 on February 7 2004 at 7:46 UT
            JD = 2453042.82222;
            for (i = 0; i < 10; i++)
            {
                GalileanMoonsDetails GalileanDetails1 = GalileanMoons.Calculate(JD);
                JD += (1.0 / 1440);
            }

            SaturnRingDetails saturnrings = SaturnRings.Calculate(2448972.50068);

            SaturnMoonsDetails saturnMoons = SaturnMoons.Calculate(2451439.50074);


            PhysicalMoonDetails MoonDetails = PhysicalMoon.CalculateGeocentric(2448724.5);
            PhysicalMoonDetails MoonDetail2 = PhysicalMoon.CalculateTopocentric(2448724.5, 10, 52);
            SelenographicMoonDetails CAASelenographicMoonDetails =
                PhysicalMoon.CalculateSelenographicPositionOfSun(2448724.5);

            double AltitudeOfSun = PhysicalMoon.AltitudeOfSun(2448724.5, -20, 9.7);
            double TimeOfSunrise = PhysicalMoon.TimeOfSunrise(2448724.5, -20, 9.7);
            double TimeOfSunset = PhysicalMoon.TimeOfSunset(2448724.5, -20, 9.7);
        }

        [Test]
        public void IlluminatedFractionTest()
        {
            double k = IlluminatedFraction.ConvertToIlluminatedFraction(0.724604, 0.983824, 0.910947);
            double pa1 = IlluminatedFraction.PhaseAngle(0.724604, 0.983824, 0.910947);
            double pa = IlluminatedFraction.PhaseAngle(0.724604, 0.983824, -2.62070, 26.11428, 88.35704, 0.910947);
            double k2 = IlluminatedFraction.ConvertToIlluminatedFraction(pa);
            double pa2 = IlluminatedFraction.PhaseAngleRectangular(0.621746, -0.664810, -0.033134, -2.62070,
                                                                   26.11428, 0.910947);
            double k3 = IlluminatedFraction.ConvertToIlluminatedFraction(pa2);

            double VenusMag = IlluminatedFraction.VenusMagnitudeMuller(0.724604, 0.910947, 72.96);
            double VenusMag2 = IlluminatedFraction.VenusMagnitudeAA(0.724604, 0.910947, 72.96);

            double SaturnMag = IlluminatedFraction.SaturnMagnitudeMuller(9.867882, 10.464606, 4.198, 16.442);
            double SaturnMag2 = IlluminatedFraction.SaturnMagnitudeAA(9.867882, 10.464606, 4.198, 16.442);


        }

        [Test]
        public void ParallaxTest()
        {
            Coordinate2D TopocentricDelta =
                Parallax.Equatorial2TopocentricDelta(CoordinateTransformation.DMSToDegrees(22, 38, 7.25), -15.771083,
                                                     0.37276, CoordinateTransformation.DMSToDegrees(7, 47, 27) * 15,
                                                     CoordinateTransformation.DMSToDegrees(33, 21, 22), 1706,
                                                     2452879.63681);
            Coordinate2D Topocentric =
                Parallax.Equatorial2Topocentric(CoordinateTransformation.DMSToDegrees(22, 38, 7.25), -15.771083,
                                                0.37276, CoordinateTransformation.DMSToDegrees(7, 47, 27) * 15,
                                                CoordinateTransformation.DMSToDegrees(33, 21, 22), 1706,
                                                2452879.63681);


            double distance2 = Parallax.ParallaxToDistance(CoordinateTransformation.DMSToDegrees(0, 59, 27.7));
            double parallax2 = Parallax.DistanceToParallax(distance2);

            TopocentricEclipticDetails TopocentricDetails =
                Parallax.Ecliptic2Topocentric(CoordinateTransformation.DMSToDegrees(181, 46, 22.5),
                                              CoordinateTransformation.DMSToDegrees(2, 17, 26.2),
                                              CoordinateTransformation.DMSToDegrees(0, 16, 15.5),
                                              Parallax.ParallaxToDistance(CoordinateTransformation.DMSToDegrees(0,
                                                                                                                59,
                                                                                                                27.7)),
                                              CoordinateTransformation.DMSToDegrees(23, 28, 0.8),
                                              0, CoordinateTransformation.DMSToDegrees(50, 5, 7.8), 0,
                                              2452879.150858);
        }

        [Test]
        public void PlanetaryPhenomenaTest()
        {
            double Kpp = AAPlus.PlanetaryPhenomena.K(1993.75, PlanetaryObject.MERCURY,
                                              EventType.INFERIOR_CONJUNCTION);
            double MercuryInferiorConjunction = AAPlus.PlanetaryPhenomena.Mean(Kpp, PlanetaryObject.MERCURY,
                                                                        EventType.INFERIOR_CONJUNCTION);
            double MercuryInferiorConjunction2 = AAPlus.PlanetaryPhenomena.True(Kpp, PlanetaryObject.MERCURY,
                                                                         EventType.INFERIOR_CONJUNCTION);

            double Kpp2 = AAPlus.PlanetaryPhenomena.K(2125.5, PlanetaryObject.SATURN, EventType.CONJUNCTION);
            double SaturnConjunction = AAPlus.PlanetaryPhenomena.Mean(Kpp2, PlanetaryObject.SATURN,
                                                               EventType.CONJUNCTION);
            double SaturnConjunction2 = AAPlus.PlanetaryPhenomena.True(Kpp2, PlanetaryObject.SATURN,
                                                                EventType.CONJUNCTION);

            double MercuryWesternElongation = AAPlus.PlanetaryPhenomena.True(Kpp, PlanetaryObject.MERCURY,
                                                                      EventType.WESTERN_ELONGATION);
            double MercuryWesternElongationValue = AAPlus.PlanetaryPhenomena.ElongationValue(Kpp,
                                                                                      PlanetaryObject.MERCURY,
                                                                                      false);

            double MarsStation2 = AAPlus.PlanetaryPhenomena.True(-2, PlanetaryObject.MARS, EventType.STATION2);

            double MercuryK = AAPlus.PlanetaryPhenomena.K(1631.8, PlanetaryObject.MERCURY,
                                                   EventType.INFERIOR_CONJUNCTION);
            double MercuryIC = AAPlus.PlanetaryPhenomena.True(MercuryK, PlanetaryObject.MERCURY,
                                                       EventType.INFERIOR_CONJUNCTION);

            double VenusKpp = AAPlus.PlanetaryPhenomena.K(1882.9, PlanetaryObject.VENUS,
                                                   EventType.INFERIOR_CONJUNCTION);
            double VenusIC = AAPlus.PlanetaryPhenomena.True(VenusKpp, PlanetaryObject.VENUS,
                                                     EventType.INFERIOR_CONJUNCTION);

            double MarsKpp = AAPlus.PlanetaryPhenomena.K(2729.65, PlanetaryObject.MARS, EventType.OPPOSITION);
            double MarsOP = AAPlus.PlanetaryPhenomena.True(MarsKpp, PlanetaryObject.MARS, EventType.OPPOSITION);

            double JupiterKpp = AAPlus.PlanetaryPhenomena.K(-5, PlanetaryObject.JUPITER, EventType.OPPOSITION);
            double JupiterOP = AAPlus.PlanetaryPhenomena.True(JupiterKpp, PlanetaryObject.JUPITER,
                                                       EventType.OPPOSITION);

            double SaturnKpp = AAPlus.PlanetaryPhenomena.K(-5, PlanetaryObject.SATURN, EventType.OPPOSITION);
            double SaturnOP = AAPlus.PlanetaryPhenomena.True(SaturnKpp, PlanetaryObject.SATURN,
                                                      EventType.OPPOSITION);

            double UranusKpp = AAPlus.PlanetaryPhenomena.K(1780.6, PlanetaryObject.URANUS, EventType.OPPOSITION);
            double UranusOP = AAPlus.PlanetaryPhenomena.True(UranusKpp, PlanetaryObject.URANUS,
                                                      EventType.OPPOSITION);

            double NeptuneKpp = AAPlus.PlanetaryPhenomena.K(1846.5, PlanetaryObject.NEPTUNE,
                                                     EventType.OPPOSITION);
            double NeptuneOP = AAPlus.PlanetaryPhenomena.True(NeptuneKpp, PlanetaryObject.NEPTUNE,
                                                       EventType.OPPOSITION);
        }

        [Test]
        public void DiametersTest()
        {
            double sd1 = Diameters.SunSemidiameterA(1);
            double sd2 = Diameters.SunSemidiameterA(2);

            double sd3 = Diameters.VenusSemidiameterA(1);
            double sd4 = Diameters.VenusSemidiameterA(2);
            double sd5 = Diameters.VenusSemidiameterB(1);
            double sd6 = Diameters.VenusSemidiameterB(2);

            double sd11 = Diameters.MarsSemidiameterA(1);
            double sd12 = Diameters.MarsSemidiameterA(2);
            double sd13 = Diameters.MarsSemidiameterB(1);
            double sd14 = Diameters.MarsSemidiameterB(2);

            double sd15 = Diameters.JupiterEquatorialSemidiameterA(1);
            double sd16 = Diameters.JupiterEquatorialSemidiameterA(2);
            double sd17 = Diameters.JupiterEquatorialSemidiameterB(1);
            double sd18 = Diameters.JupiterEquatorialSemidiameterB(2);

            double sd19 = Diameters.JupiterPolarSemidiameterA(1);
            double sd20 = Diameters.JupiterPolarSemidiameterA(2);
            double sd21 = Diameters.JupiterPolarSemidiameterB(1);
            double sd22 = Diameters.JupiterPolarSemidiameterB(2);

            double sd23 = Diameters.SaturnEquatorialSemidiameterA(1);
            double sd24 = Diameters.SaturnEquatorialSemidiameterA(2);
            double sd25 = Diameters.SaturnEquatorialSemidiameterB(1);
            double sd26 = Diameters.SaturnEquatorialSemidiameterB(2);

            double sd27 = Diameters.SaturnPolarSemidiameterA(1);
            double sd28 = Diameters.SaturnPolarSemidiameterA(2);
            double sd29 = Diameters.SaturnPolarSemidiameterB(1);
            double sd30 = Diameters.SaturnPolarSemidiameterB(2);

            double sd31 = Diameters.ApparentSaturnPolarSemidiameterA(1, 16.442);
            double sd32 = Diameters.ApparentSaturnPolarSemidiameterA(2, 16.442);

            double sd33 = Diameters.UranusSemidiameterA(1);
            double sd34 = Diameters.UranusSemidiameterA(2);
            double sd35 = Diameters.UranusSemidiameterB(1);
            double sd36 = Diameters.UranusSemidiameterB(2);

            double sd37 = Diameters.NeptuneSemidiameterA(1);
            double sd38 = Diameters.NeptuneSemidiameterA(2);
            double sd39 = Diameters.NeptuneSemidiameterB(1);
            double sd40 = Diameters.NeptuneSemidiameterB(2);

            double sd41 = Diameters.PlutoSemidiameterB(1);
            double sd42 = Diameters.PlutoSemidiameterB(2);

            double sd43 = Diameters.GeocentricMoonSemidiameter(368407.9);
            double sd44 = Diameters.GeocentricMoonSemidiameter(368407.9 - 10000);

            double sd45 = Diameters.TopocentricMoonSemidiameter(368407.9, 5, 0, 33.356111, 1706);
            double sd46 = Diameters.TopocentricMoonSemidiameter(368407.9, 5, 6, 33.356111, 1706);
            double sd47 = Diameters.TopocentricMoonSemidiameter(368407.9 - 10000, 5, 0, 33.356111, 1706);
            double sd48 = Diameters.TopocentricMoonSemidiameter(368407.9 - 10000, 5, 6, 33.356111, 1706);

            double sd49 = Diameters.AsteroidDiameter(4, 0.04);
            double sd50 = Diameters.AsteroidDiameter(4, 0.08);
            double sd51 = Diameters.AsteroidDiameter(6, 0.04);
            double sd53 = Diameters.AsteroidDiameter(6, 0.08);
            double sd54 = Diameters.ApparentAsteroidDiameter(1, 250);
            double sd55 = Diameters.ApparentAsteroidDiameter(1, 1000);
        }

        [Test]
        public void RiseTransitSetTest()
        {
            long Year = 0;
            long Month = 0;
            long Day = 0;
            double Alpha1 = CoordinateTransformation.DMSToDegrees(2, 42, 43.25);
            double Alpha2 = CoordinateTransformation.DMSToDegrees(2, 46, 55.51);
            double Alpha3 = CoordinateTransformation.DMSToDegrees(2, 51, 07.69);
            double Delta1 = CoordinateTransformation.DMSToDegrees(18, 02, 51.4);
            double Delta2 = CoordinateTransformation.DMSToDegrees(18, 26, 27.3);
            double Delta3 = CoordinateTransformation.DMSToDegrees(18, 49, 38.7);
            double JD2 = 2447240.5;
            double Longitude = 71.0833;
            double Latitude = 42.3333;
            RiseTransitSetDetails RiseTransitSetTime = RiseTransitSet.Rise(JD2, Alpha1, Delta1, Alpha2, Delta2,
                                                                           Alpha3, Delta3, Longitude, Latitude,
                                                                           -0.5667);
            {
                if (RiseTransitSetTime.bValid)
                {
                    double riseJD = (JD2 + (RiseTransitSetTime.Rise / 24.00));
                    Date rtsDate = new Date(riseJD, true);
                    long Hours;
                    long Minutes;
                    double Sec;
                    rtsDate.Get(out Year, out Month, out Day, out Hours, out Minutes, out Sec);
                    Printf("Venus rise for Boston for UTC %d/%d/%d occurs at %02d:%02d:%02d\n", Year,
                           Month, Day, Hours, Minutes, (int)(Sec));
                    double transitJD = (JD2 + (RiseTransitSetTime.Transit / 24.00));
                    rtsDate = new Date(transitJD, true);
                    rtsDate.Get(out Year, out Month, out Day, out Hours, out Minutes, out Sec);
                    Printf("Venus transit for Boston for UTC %d/%d/%d occurs at %02d:%02d:%02d\n", Year,
                           Month, Day, Hours, Minutes, (int)(Sec));
                    double setJD = (JD2 + (RiseTransitSetTime.Set / 24.00));
                    rtsDate = new Date(setJD, true);
                    rtsDate.Get(out Year, out Month, out Day, out Hours, out Minutes, out Sec);
                    Printf("Venus set for Boston UTC %d/%d/%d occurs at %02d:%02d:%02d\n", Year, Month,
                           Day,
                           Hours, Minutes, (int)(Sec));
                }
                else
                {
                    Printf("Venus does not rise, transit or set for Boston for UTC %d/%d/%d\n", Year,
                           Month,
                           Day);
                }
            }
            {

                //Calculate the time of moon set for 11th of August 2009 UTC for Palamor Observatory 
                int YYYY = 2009;
                int MM = 8;
                int DD = 11;
                Date CalcDate = new Date(YYYY, MM, DD, true);

                JD2 = CalcDate.Julian();
                LunarRaDec(JD2 - 1, out Alpha1, out Delta1);
                LunarRaDec(JD2, out Alpha2, out Delta2);
                LunarRaDec(JD2 + 1, out Alpha3, out Delta3);
                Longitude = CoordinateTransformation.DMSToDegrees(116, 51, 45); //West is considered positive
                Latitude = CoordinateTransformation.DMSToDegrees(33, 21, 22);
                RiseTransitSetTime = RiseTransitSet.Rise(JD2, Alpha1, Delta1, Alpha2, Delta2, Alpha3, Delta3,
                                                         Longitude,
                                                         Latitude, -0.8333);
                if (RiseTransitSetTime.bValid)
                {
                    double riseJD = (JD2 + (RiseTransitSetTime.Rise / 24.00));
                    Date rtsDate = new Date(riseJD, true);
                    long Hours;
                    long Minutes;
                    double Sec;
                    rtsDate.Get(out Year, out Month, out Day, out Hours, out Minutes, out Sec);
                    Printf(
                        "Moon rise for Palamor Observatory for UTC %d/%d/%d occurs at %02d:%02d:%02d\n",
                        Year, Month, Day, Hours, Minutes, (int)(Sec));
                    double transitJD = (JD2 + (RiseTransitSetTime.Transit / 24.00));
                    rtsDate = new Date(transitJD, true);
                    rtsDate.Get(out Year, out Month, out Day, out Hours, out Minutes, out Sec);
                    Printf(
                        "Moon transit for Palamor Observatory for UTC %d/%d/%d occurs at %02d:%02d:%02d\n", Year,
                        Month,
                        Day, Hours, Minutes, (int)(Sec));
                    double setJD = (JD2 + (RiseTransitSetTime.Set / 24.00));
                    rtsDate = new Date(setJD, true);
                    rtsDate.Get(out Year, out Month, out Day, out Hours, out Minutes, out Sec);
                    Printf(
                        "Moon set for Palamor Observatory for UTC %d/%d/%d occurs at %02d:%02d:%02d\n",
                        Year, Month, Day, Hours, Minutes, (int)(Sec));
                }
                else
                {
                    Printf(
                        "Moon for Palamor Observatory does not rise, transit or set for UTC %d/%d/%d\n",
                        Year, Month, Day);
                }
            }
        }

        [Test]
        public void ParabolicTest()
        {
            ParabolicObjectElements elements2 = new ParabolicObjectElements(); //Elements taken from http://www.cfa.harvard.edu/mpec/J98/J98H29.html
            elements2.q = 1.48678;
            elements2.i = 104.668; //J2000
            elements2.omega = 222.103; //J2000
            elements2.w = 1.146; //J2000
            elements2.T = Date.DateToJD(1998, 4, 14.205, true);
            elements2.JDEquinox = 2451545.0; //J2000
            ParabolicObjectDetails details2 = Parabolic.Calculate(2451030.5, elements2);

            EllipticalObjectElements elements3 = new EllipticalObjectElements();
            elements3.a = 17.9400782;
            elements3.e = 0.96727426;
            elements3.i = 0; //Not required
            elements3.omega = 0; //Not required
            elements3.w = 111.84644;
            elements3.T = 2446470.5 + 0.45891;
            elements3.JDEquinox = 0; //Not required
            NodeObjectDetails nodedetails = Nodes.PassageThroAscendingNode(elements3);
            NodeObjectDetails nodedetails2 = Nodes.PassageThroDescendingNode(elements3);

            ParabolicObjectElements elements4 = new ParabolicObjectElements();
            elements4.q = 1.324502;
            elements4.i = 0; //Not required
            elements4.omega = 0; //Not required
            elements4.w = 154.9103;
            elements4.T = 2447758.5 + 0.2910;
            elements4.JDEquinox = 0; //Not required
            NodeObjectDetails nodedetails3 = Nodes.PassageThroAscendingNode(elements4);
            NodeObjectDetails nodedetails4 = Nodes.PassageThroDescendingNode(elements4);


            EllipticalObjectElements elements5 = new EllipticalObjectElements();
            elements5.a = 0.723329820;
            elements5.e = 0.00678195;
            elements5.i = 0; //Not required
            elements5.omega = 0; //Not required
            elements5.w = 54.778485;
            elements5.T = 2443873.704;
            elements5.JDEquinox = 0; //Not required
            NodeObjectDetails nodedetails5 = Nodes.PassageThroAscendingNode(elements5);

            double MoonK2 = MoonNodes.K(1987.37);
            double MoonJD = MoonNodes.PassageThroNode(-170);


            double Y = Interpolation.Interpolate(0.18125, 0.884226, 0.877366, 0.870531);

            double NM = 0d;
            double YM = Interpolation.Extremum(1.3814294, 1.3812213, 1.3812453, out NM);

            double N0 = Interpolation.Zero(-1693.4, 406.3, 2303.2);

            double N02 = Interpolation.Zero2(-2, 3, 2);

            double Y2 = Interpolation.Interpolate(0.2777778, 36.125, 24.606, 15.486, 8.694, 4.133);

            double N03 = Interpolation.Zero(CoordinateTransformation.DMSToDegrees(1, 11, 21.23, false),
                                            CoordinateTransformation.DMSToDegrees(0, 28, 12.31, false),
                                            CoordinateTransformation.DMSToDegrees(0, 16, 7.02),
                                            CoordinateTransformation.DMSToDegrees(1, 1, 0.13),
                                            CoordinateTransformation.DMSToDegrees(1, 45, 46.33));

            double N04 = Interpolation.Zero(CoordinateTransformation.DMSToDegrees(0, 28, 12.31, false),
                                            CoordinateTransformation.DMSToDegrees(0, 16, 7.02),
                                            CoordinateTransformation.DMSToDegrees(1, 1, 0.13));

            double N05 = Interpolation.Zero2(-13, -2, 3, 2, -5);

            double Y3 = Interpolation.InterpolateToHalves(1128.732, 1402.835, 1677.247, 1951.983);

            double[] X1 =
                {
                    29.43,
                    30.97,
                    27.69,
                    28.11,
                    31.58,
                    33.05
                }
                ;
            double[] Y1 =
                {
                    0.4913598528,
                    0.5145891926,
                    0.4646875083,
                    0.4711658342,
                    0.5236885653,
                    0.5453707057
                };

            double Y4 = Interpolation.LagrangeInterpolate(30, 6, X1, Y1);
            double Y5 = Interpolation.LagrangeInterpolate(0, 6, X1, Y1);
            double Y6 = Interpolation.LagrangeInterpolate(90, 6, X1, Y1);
        }

        [Test]
        public void EllipticalTest()
        {
            EllipticalPlanetaryDetails VenusDetails = Elliptical.Calculate(2448976.5, Elliptical.EllipticalObject.Venus);

            EllipticalPlanetaryDetails SunDetails = Elliptical.Calculate(2453149.5, Elliptical.EllipticalObject.Sun);

            EllipticalObjectElements elements = new EllipticalObjectElements();
            elements.a = 2.2091404;
            elements.e = 0.8502196;
            elements.i = 11.94524;
            elements.omega = 334.75006;
            elements.w = 186.23352;
            elements.T = 2448192.5 + 0.54502;
            elements.JDEquinox = 2451544.5;
            EllipticalObjectDetails details = Elliptical.Calculate(2448170.5, out elements);

            double Velocity1 = Elliptical.InstantaneousVelocity(1, 17.9400782);
            double Velocity2 = Elliptical.VelocityAtPerihelion(0.96727426, 17.9400782);
            double Velocity3 = Elliptical.VelocityAtAphelion(0.96727426, 17.9400782);

            double Length = Elliptical.LengthOfEllipse(0.96727426, 17.9400782);

            double Mag1 = Elliptical.MinorPlanetMagnitude(3.34, 1.6906631928, 0.12, 2.6154983761, 120);
            double Mag2 = Elliptical.CometMagnitude(5.5, 0.378, 10, 0.658);
            double Mag3 = Elliptical.CometMagnitude(5.5, 1.1017, 10, 1.5228);

        }

        [Test]
        public void MoonIlluminatedFractionTest()
        {
            double MoonGeocentricElongation = AAPlus.MoonIlluminatedFraction.GeocentricElongation(8.97922, 13.7684,
                                                                                                  1.377194, 8.6964);
            double MoonPhaseAngle = AAPlus.MoonIlluminatedFraction.PhaseAngle(MoonGeocentricElongation, 368410, 149971520);
            double MoonIlluminatedFraction = AAPlus.MoonIlluminatedFraction.IlluminatedFraction(MoonPhaseAngle);
            double MoonPositionAngle =
                AAPlus.MoonIlluminatedFraction.PositionAngle(CoordinateTransformation.DMSToDegrees(1, 22, 37.9), 8.6964,
                                                             134.6885 / 15, 13.7684);
        }

        [Test]
        public void ElementsPlanetaryOrbitTest()
        {
            double Mer_L = ElementsPlanetaryOrbit.MercuryMeanLongitude(2475460.5);
            double Mer_a = ElementsPlanetaryOrbit.MercurySemimajorAxis(2475460.5);
            double Mer_e = ElementsPlanetaryOrbit.MercuryEccentricity(2475460.5);
            double Mer_i = ElementsPlanetaryOrbit.MercuryInclination(2475460.5);
            double Mer_omega = ElementsPlanetaryOrbit.MercuryLongitudeAscendingNode(2475460.5);
            double mer_pi = ElementsPlanetaryOrbit.MercuryLongitudePerihelion(2475460.5);

            double Ven_L = ElementsPlanetaryOrbit.VenusMeanLongitude(2475460.5);
            double Ven_a = ElementsPlanetaryOrbit.VenusSemimajorAxis(2475460.5);
            double Ven_e = ElementsPlanetaryOrbit.VenusEccentricity(2475460.5);
            double Ven_i = ElementsPlanetaryOrbit.VenusInclination(2475460.5);
            double Ven_omega = ElementsPlanetaryOrbit.VenusLongitudeAscendingNode(2475460.5);
            double Ven_pi = ElementsPlanetaryOrbit.VenusLongitudePerihelion(2475460.5);

            double Ea_L = ElementsPlanetaryOrbit.EarthMeanLongitude(2475460.5);
            double Ea_a = ElementsPlanetaryOrbit.EarthSemimajorAxis(2475460.5);
            double Ea_e = ElementsPlanetaryOrbit.EarthEccentricity(2475460.5);
            double Ea_i = ElementsPlanetaryOrbit.EarthInclination(2475460.5);
            double Ea_pi = ElementsPlanetaryOrbit.EarthLongitudePerihelion(2475460.5);

            double Mars_L = ElementsPlanetaryOrbit.MarsMeanLongitude(2475460.5);
            double Mars_a = ElementsPlanetaryOrbit.MarsSemimajorAxis(2475460.5);
            double Mars_e = ElementsPlanetaryOrbit.MarsEccentricity(2475460.5);
            double Mars_i = ElementsPlanetaryOrbit.MarsInclination(2475460.5);
            double Mars_omega = ElementsPlanetaryOrbit.MarsLongitudeAscendingNode(2475460.5);
            double Mars_pi = ElementsPlanetaryOrbit.MarsLongitudePerihelion(2475460.5);

            double Jup_L = ElementsPlanetaryOrbit.JupiterMeanLongitude(2475460.5);
            double Jup_a = ElementsPlanetaryOrbit.JupiterSemimajorAxis(2475460.5);
            double Jup_e = ElementsPlanetaryOrbit.JupiterEccentricity(2475460.5);
            double Jup_i = ElementsPlanetaryOrbit.JupiterInclination(2475460.5);
            double Jup_omega = ElementsPlanetaryOrbit.JupiterLongitudeAscendingNode(2475460.5);
            double Jup_pi = ElementsPlanetaryOrbit.JupiterLongitudePerihelion(2475460.5);

            double Sat_L = ElementsPlanetaryOrbit.SaturnMeanLongitude(2475460.5);
            double Sat_a = ElementsPlanetaryOrbit.SaturnSemimajorAxis(2475460.5);
            double Sat_e = ElementsPlanetaryOrbit.SaturnEccentricity(2475460.5);
            double Sat_i = ElementsPlanetaryOrbit.SaturnInclination(2475460.5);
            double Sat_omega = ElementsPlanetaryOrbit.SaturnLongitudeAscendingNode(2475460.5);
            double Sat_pi = ElementsPlanetaryOrbit.SaturnLongitudePerihelion(2475460.5);

            double Ura_L = ElementsPlanetaryOrbit.UranusMeanLongitude(2475460.5);
            double Ura_a = ElementsPlanetaryOrbit.UranusSemimajorAxis(2475460.5);
            double Ura_e = ElementsPlanetaryOrbit.UranusEccentricity(2475460.5);
            double Ura_i = ElementsPlanetaryOrbit.UranusInclination(2475460.5);
            double Ura_omega = ElementsPlanetaryOrbit.UranusLongitudeAscendingNode(2475460.5);
            double Ura_pi = ElementsPlanetaryOrbit.UranusLongitudePerihelion(2475460.5);

            double Nep_L = ElementsPlanetaryOrbit.NeptuneMeanLongitude(2475460.5);
            double Nep_a = ElementsPlanetaryOrbit.NeptuneSemimajorAxis(2475460.5);
            double Nep_e = ElementsPlanetaryOrbit.NeptuneEccentricity(2475460.5);
            double Nep_i = ElementsPlanetaryOrbit.NeptuneInclination(2475460.5);
            double Nep_omega = ElementsPlanetaryOrbit.NeptuneLongitudeAscendingNode(2475460.5);
            double Nep_pi = ElementsPlanetaryOrbit.NeptuneLongitudePerihelion(2475460.5);


            double Mer_L2 = ElementsPlanetaryOrbit.MercuryMeanLongitudeJ2000(2475460.5);
            double Mer_i2 = ElementsPlanetaryOrbit.MercuryInclinationJ2000(2475460.5);
            double Mer_omega2 = ElementsPlanetaryOrbit.MercuryLongitudeAscendingNodeJ2000(2475460.5);
            double mer_pi2 = ElementsPlanetaryOrbit.MercuryLongitudePerihelionJ2000(2475460.5);

            double Ven_L2 = ElementsPlanetaryOrbit.VenusMeanLongitudeJ2000(2475460.5);
            double Ven_i2 = ElementsPlanetaryOrbit.VenusInclinationJ2000(2475460.5);
            double Ven_omega2 = ElementsPlanetaryOrbit.VenusLongitudeAscendingNodeJ2000(2475460.5);
            double Ven_pi2 = ElementsPlanetaryOrbit.VenusLongitudePerihelionJ2000(2475460.5);

            double Ea_L2 = ElementsPlanetaryOrbit.EarthMeanLongitudeJ2000(2475460.5);
            double Ea_i2 = ElementsPlanetaryOrbit.EarthInclinationJ2000(2475460.5);
            double Ea_omega2 = ElementsPlanetaryOrbit.EarthLongitudeAscendingNodeJ2000(2475460.5);
            double Ea_pi2 = ElementsPlanetaryOrbit.EarthLongitudePerihelionJ2000(2475460.5);

            double Mars_L2 = ElementsPlanetaryOrbit.MarsMeanLongitudeJ2000(2475460.5);
            double Mars_i2 = ElementsPlanetaryOrbit.MarsInclinationJ2000(2475460.5);
            double Mars_omega2 = ElementsPlanetaryOrbit.MarsLongitudeAscendingNodeJ2000(2475460.5);
            double Mars_pi2 = ElementsPlanetaryOrbit.MarsLongitudePerihelionJ2000(2475460.5);

            double Jup_L2 = ElementsPlanetaryOrbit.JupiterMeanLongitudeJ2000(2475460.5);
            double Jup_i2 = ElementsPlanetaryOrbit.JupiterInclinationJ2000(2475460.5);
            double Jup_omega2 = ElementsPlanetaryOrbit.JupiterLongitudeAscendingNodeJ2000(2475460.5);
            double Jup_pi2 = ElementsPlanetaryOrbit.JupiterLongitudePerihelionJ2000(2475460.5);

            double Sat_L2 = ElementsPlanetaryOrbit.SaturnMeanLongitudeJ2000(2475460.5);
            double Sat_i2 = ElementsPlanetaryOrbit.SaturnInclinationJ2000(2475460.5);
            double Sat_omega2 = ElementsPlanetaryOrbit.SaturnLongitudeAscendingNodeJ2000(2475460.5);
            double Sat_pi2 = ElementsPlanetaryOrbit.SaturnLongitudePerihelionJ2000(2475460.5);

            double Ura_L2 = ElementsPlanetaryOrbit.UranusMeanLongitudeJ2000(2475460.5);
            double Ura_i2 = ElementsPlanetaryOrbit.UranusInclinationJ2000(2475460.5);
            double Ura_omega2 = ElementsPlanetaryOrbit.UranusLongitudeAscendingNodeJ2000(2475460.5);
            double Ura_pi2 = ElementsPlanetaryOrbit.UranusLongitudePerihelionJ2000(2475460.5);

            double Nep_L2 = ElementsPlanetaryOrbit.NeptuneMeanLongitudeJ2000(2475460.5);
            double Nep_i2 = ElementsPlanetaryOrbit.NeptuneInclinationJ2000(2475460.5);
            double Nep_omega2 = ElementsPlanetaryOrbit.NeptuneLongitudeAscendingNodeJ2000(2475460.5);
            double Nep_pi2 = ElementsPlanetaryOrbit.NeptuneLongitudePerihelionJ2000(2475460.5);
        }

        [Test]
        public void EquationOfTimeTest()
        {
            double E = EquationOfTime.Calculate(2448908.5);
        }

        [Test]
        public void EclipticalElementsTest()
        {
            EclipticalElementDetails ed1 = AAPlus.EclipticalElements.Calculate(47.1220, 151.4486, 45.7481, 2358042.5305,
                                                                        2433282.4235);
            EclipticalElementDetails ed2 = AAPlus.EclipticalElements.Calculate(11.93911, 186.24444, 334.04096, 2433282.4235,
                                                                        2451545.0);
            EclipticalElementDetails ed3 = AAPlus.EclipticalElements.FK4B1950ToFK5J2000(11.93911, 186.24444, 334.04096);
            EclipticalElementDetails ed4 = AAPlus.EclipticalElements.FK4B1950ToFK5J2000(145, 186.24444, 334.04096);
        }

        [Test]
        public void MoonPerigeeApogeeTest()
        {
            double MoonK = AAPlus.MoonPerigeeApogee.K(1988.75);
            double MoonApogee = AAPlus.MoonPerigeeApogee.MeanApogee(-148.5);
            double MoonApogee2 = AAPlus.MoonPerigeeApogee.TrueApogee(-148.5);
            double MoonApogeeParallax = AAPlus.MoonPerigeeApogee.ApogeeParallax(-148.5);
            double MoonApogeeDistance = Moon.HorizontalParallaxToRadiusVector(MoonApogeeParallax);

            MoonK = AAPlus.MoonPerigeeApogee.K(1990.9);
            double MoonPerigee = AAPlus.MoonPerigeeApogee.MeanPerigee(-120);
            double MoonPerigee2 = AAPlus.MoonPerigeeApogee.TruePerigee(-120);
            MoonK = AAPlus.MoonPerigeeApogee.K(1930.0);
            double MoonPerigee3 = AAPlus.MoonPerigeeApogee.TruePerigee(-927);
            double MoonPerigeeParallax = AAPlus.MoonPerigeeApogee.PerigeeParallax(-927);
            double MoonRadiusVector = Moon.HorizontalParallaxToRadiusVector(MoonPerigeeParallax);
            double MoonRadiusVector2 = Moon.HorizontalParallaxToRadiusVector(0.991990);
            double MoonParallax2 = Moon.RadiusVectorToHorizontalParallax(MoonRadiusVector2);
        }

        [Test]
        public void PlanetPerihelionAphelionTest()
        {
            long VenusK = AAPlus.PlanetPerihelionAphelion.VenusK(1978.79);
            double VenusPerihelion = AAPlus.PlanetPerihelionAphelion.VenusPerihelion(VenusK);

            long MarsK = AAPlus.PlanetPerihelionAphelion.MarsK(2032);
            double MarsAphelion = AAPlus.PlanetPerihelionAphelion.MarsAphelion(MarsK);

            long SaturnK = AAPlus.PlanetPerihelionAphelion.SaturnK(1925);
            double SaturnAphelion = AAPlus.PlanetPerihelionAphelion.SaturnAphelion(SaturnK);
            SaturnK = AAPlus.PlanetPerihelionAphelion.SaturnK(1940);
            double SaturnPerihelion = AAPlus.PlanetPerihelionAphelion.SaturnPerihelion(SaturnK);

            long UranusK = AAPlus.PlanetPerihelionAphelion.UranusK(1750);
            double UranusAphelion = AAPlus.PlanetPerihelionAphelion.UranusAphelion(UranusK);
            UranusK = AAPlus.PlanetPerihelionAphelion.UranusK(1890);
            double UranusPerihelion = AAPlus.PlanetPerihelionAphelion.UranusPerihelion(UranusK);
            UranusK = AAPlus.PlanetPerihelionAphelion.UranusK(2060);
            UranusPerihelion = AAPlus.PlanetPerihelionAphelion.UranusPerihelion(UranusK);

            double EarthPerihelion = AAPlus.PlanetPerihelionAphelion.EarthPerihelion(-10, true);
            double EarthPerihelion2 = AAPlus.PlanetPerihelionAphelion.EarthPerihelion(-10, false);
        }

        [Test]
        public void MoonTest()
        {
            double MoonLong;
            double MoonLat;
            MoonLong = Moon.EclipticLongitude(2448724.5);
            MoonLat = Moon.EclipticLatitude(2448724.5);
            double MoonRadius = Moon.RadiusVector(2448724.5);
            double MoonParallax = Moon.RadiusVectorToHorizontalParallax(MoonRadius);
            double MoonMeanAscendingNode = Moon.MeanLongitudeAscendingNode(2448724.5);
            double TrueMeanAscendingNode = Moon.TrueLongitudeAscendingNode(2448724.5);
            double MoonMeanPerigee = Moon.MeanLongitudePerigee(2448724.5);
        }

        [Test]
        public void PlutoTest()
        {
            double PlutoLong = Pluto.EclipticLongitude(2448908.5);
            double PlutoLat = Pluto.EclipticLatitude(2448908.5);
            double PlutoRadius = Pluto.RadiusVector(2448908.5);
        }

        [Test]
        public void SaturnTest()
        {
            double SaturnLong = Saturn.EclipticLongitude(2448972.50068);
            double SaturnLat = Saturn.EclipticLatitude(2448972.50068);
            double SaturnRadius = Saturn.RadiusVector(2448972.50068);
        }

        [Test]
        public void UranusTest()
        {
            double UranusLong = Uranus.EclipticLongitude(2448976.5);
            double UranusLat = Uranus.EclipticLatitude(2448976.5);
            double UranusRadius = Uranus.RadiusVector(2448976.5);
        }

        [Test]
        public void NeptuneTest()
        {
            double NeptuneLong = Neptune.EclipticLongitude(2448935.500683);
            double NeptuneLat = Neptune.EclipticLatitude(2448935.500683);
            double NeptuneRadius = Neptune.RadiusVector(2448935.500683);
        }

        [Test]
        public void JupiterTest()
        {
            double JupiterLong = Jupiter.EclipticLongitude(2448972.50068);
            double JupiterLat = Jupiter.EclipticLatitude(2448972.50068);
            double JupiterRadius = Jupiter.RadiusVector(2448972.50068);
        }

        [Test]
        public void MarsTest()
        {
            double MarsLong = Mars.EclipticLongitude(2448935.500683);
            double MarsLat = Mars.EclipticLatitude(2448935.500683);
            double MarsRadius = Mars.RadiusVector(2448935.500683);
        }

        [Test]
        public void EquinoxesAndSolstices()
        {
            long Month;
            long Day;
            Date date = new Date();
            long Year;
            long Hour;
            long Minute;
            double Second;
            double JuneSolstice = AAPlus.EquinoxesAndSolstices.SummerSolstice(1962);


            double MarchEquinox2 = AAPlus.EquinoxesAndSolstices.SpringEquinox(1996);
            date.Set(MarchEquinox2, true);
            date.Get(out Year, out Month, out Day, out Hour, out Minute, out Second);
            double JuneSolstice2 = AAPlus.EquinoxesAndSolstices.SummerSolstice(1996);
            date.Set(JuneSolstice2, true);
            date.Get(out Year, out Month, out Day, out Hour, out Minute, out Second);
            double SeptemberEquinox2 = AAPlus.EquinoxesAndSolstices.AutumnEquinox(1996);
            date.Set(SeptemberEquinox2, true);
            date.Get(out Year, out Month, out Day, out Hour, out Minute, out Second);
            double DecemberSolstice2 = AAPlus.EquinoxesAndSolstices.WinterSolstice(1996);
            date.Set(DecemberSolstice2, true);
            date.Get(out Year, out Month, out Day, out Hour, out Minute, out Second);

            DecemberSolstice2 = AAPlus.EquinoxesAndSolstices.WinterSolstice(2000);
            date.Set(DecemberSolstice2, true);
            date.Get(out Year, out Month, out Day, out Hour, out Minute, out Second);

            DecemberSolstice2 = AAPlus.EquinoxesAndSolstices.WinterSolstice(1997);
            date.Set(DecemberSolstice2, true);
            date.Get(out Year, out Month, out Day, out Hour, out Minute, out Second);

            DecemberSolstice2 = AAPlus.EquinoxesAndSolstices.WinterSolstice(2003);
            date.Set(DecemberSolstice2, true);
            date.Get(out Year, out Month, out Day, out Hour, out Minute, out Second);

            JuneSolstice2 = AAPlus.EquinoxesAndSolstices.SummerSolstice(2003);
            date.Set(JuneSolstice2, true);
            date.Get(out Year, out Month, out Day, out Hour, out Minute, out Second);

            double SpringLength = AAPlus.EquinoxesAndSolstices.LengthOfSpring(2000);
            double SummerLength = AAPlus.EquinoxesAndSolstices.LengthOfSummer(2000);
            double AutumnLength = AAPlus.EquinoxesAndSolstices.LengthOfAutumn(2000);
            double WinterLength = AAPlus.EquinoxesAndSolstices.LengthOfWinter(2000);

            SpringLength = AAPlus.EquinoxesAndSolstices.LengthOfSpring(-2000);
            SummerLength = AAPlus.EquinoxesAndSolstices.LengthOfSummer(-2000);
            AutumnLength = AAPlus.EquinoxesAndSolstices.LengthOfAutumn(-2000);
            WinterLength = AAPlus.EquinoxesAndSolstices.LengthOfWinter(-2000);

            SpringLength = AAPlus.EquinoxesAndSolstices.LengthOfSpring(4000);
            SummerLength = AAPlus.EquinoxesAndSolstices.LengthOfSummer(4000);
            AutumnLength = AAPlus.EquinoxesAndSolstices.LengthOfAutumn(4000);
            WinterLength = AAPlus.EquinoxesAndSolstices.LengthOfWinter(4000);
        }

        [Test]
        public void SunTest()
        {
            double SunLong;
            double SunLat;
            Coordinate2D Equatorial;
            SunLong = Sun.GeometricEclipticLongitude(2448908.5);
            SunLat = Sun.GeometricEclipticLatitude(2448908.5);

            double SunLongCorrection = FK5.CorrectionInLongitude(SunLong, SunLat, 2448908.5);
            double SunLatCorrection = FK5.CorrectionInLatitude(SunLong, 2448908.5);

            SunLong = Sun.ApparentEclipticLongitude(2448908.5);
            SunLat = Sun.ApparentEclipticLatitude(2448908.5);
            Equatorial = CoordinateTransformation.Ecliptic2Equatorial(SunLong, SunLat,
                                                                      Nutation.TrueObliquityOfEcliptic(2448908.5));

            Coordinate3D SunCoord = Sun.EquatorialRectangularCoordinatesMeanEquinox(2448908.5);
            Coordinate3D SunCoord2 = Sun.EclipticRectangularCoordinatesJ2000(2448908.5);
            Coordinate3D SunCoord3 = Sun.EquatorialRectangularCoordinatesJ2000(2448908.5);
            Coordinate3D SunCoord4 = Sun.EquatorialRectangularCoordinatesB1950(2448908.5);
            Coordinate3D SunCoord5 = Sun.EquatorialRectangularCoordinatesAnyEquinox(2448908.5, 2467616.0);
        }

        [Test]
        public void EarthTest()
        {
            double EarthLong = Earth.EclipticLongitude(2448908.5);
            double EarthLat = Earth.EclipticLatitude(2448908.5);
            double EarthRadius = Earth.RadiusVector(2448908.5);

            double EarthLong2 = Earth.EclipticLongitudeJ2000(2448908.5);
            double EarthLat2 = Earth.EclipticLatitudeJ2000(2448908.5);
        }

        [Test]
        public void MercuryTest()
        {
            double MercuryLong = Mercury.EclipticLongitude(2448976.5);
            double MercuryLat = Mercury.EclipticLatitude(2448976.5);
            double MercuryRadius = Mercury.RadiusVector(2448976.5);
        }

        public void VenusTest()
        {
            double VenusLong = Venus.EclipticLongitude(2448976.5);
            double VenusLat = Venus.EclipticLatitude(2448976.5);
            double VenusRadius = Venus.RadiusVector(2448976.5);
        }

        [Test]
        public void StellarMagnitudesTest()
        {
            double CombinedMag = StellarMagnitudes.CombinedMagnitude(1.96, 2.89);
            double[] Mags =
                {
                    4.73,
                    5.22,
                    5.60
                }
                ;
            double CombinedMag2 = StellarMagnitudes.CombinedMagnitude(3, Mags);
            double BrightnessRatio = StellarMagnitudes.BrightnessRatio(0.14, 2.12);
            double MagDiff = StellarMagnitudes.MagnitudeDifference(BrightnessRatio);
            double MagDiff2 = StellarMagnitudes.MagnitudeDifference(500);
        }

        [Test]
        public void BinaryStarTest()
        {
            BinaryStarDetails bsdetails = BinaryStar.Calculate(1980, 41.623, 1934.008, 0.2763, 0.907, 59.025,
                                                               23.717, 219.907);
            double ApparentE = BinaryStar.ApparentEccentricity(0.2763, 59.025, 219.907);
        }

        [Test]
        public void KeplerTest()
        {
            double E0 = Kepler.Calculate(5, 0.1, 100);
            double E02 = Kepler.Calculate(5, 0.9, 100);
            //double E03 = Kepler.Calculate(

        }

        [Test]
        public void AngularSeparationTest()
        {
            double Alpha;
            double Obliquity;
            double NutationInLongitude;
            double NutationInEcliptic;
            double AngularSeparation =
                AAPlus.AngularSeparation.Separation(CoordinateTransformation.DMSToDegrees(14, 15, 39.7),
                                                    CoordinateTransformation.DMSToDegrees(19, 10, 57),
                                                    CoordinateTransformation.DMSToDegrees(13, 25, 11.6),
                                                    CoordinateTransformation.DMSToDegrees(11, 9, 41, false));
            double AngularSeparation2 = AAPlus.AngularSeparation.Separation(
                CoordinateTransformation.DMSToDegrees(2, 0, 0), CoordinateTransformation.DMSToDegrees(0, 0, 0),
                CoordinateTransformation.DMSToDegrees(2, 0, 0), CoordinateTransformation.DMSToDegrees(0, 0, 0));
            double AngularSeparation3 = AAPlus.AngularSeparation.Separation(
                CoordinateTransformation.DMSToDegrees(2, 0, 0), CoordinateTransformation.DMSToDegrees(0, 0, 0),
                CoordinateTransformation.DMSToDegrees(14, 0, 0), CoordinateTransformation.DMSToDegrees(0, 0, 0));

            double PA0 = AAPlus.AngularSeparation.PositionAngle(CoordinateTransformation.DMSToDegrees(5, 32, 0.4),
                                                                CoordinateTransformation.DMSToDegrees(0, 17, 56.9, false),
                                                                CoordinateTransformation.DMSToDegrees(5, 36, 12.81),
                                                                CoordinateTransformation.DMSToDegrees(1, 12, 7, false));

            double PA1 = AAPlus.AngularSeparation.PositionAngle(CoordinateTransformation.DMSToDegrees(5, 40, 45.52),
                                                                CoordinateTransformation.DMSToDegrees(1, 56, 33.3, false),
                                                                CoordinateTransformation.DMSToDegrees(5, 36, 12.81),
                                                                CoordinateTransformation.DMSToDegrees(1, 12, 7, false));


            double distance =
                AAPlus.AngularSeparation.DistanceFromGreatArc(CoordinateTransformation.DMSToDegrees(5, 32, 0.4),
                                                              CoordinateTransformation.DMSToDegrees(0, 17, 56.9, false),
                                                              CoordinateTransformation.DMSToDegrees(5, 40, 45.52),
                                                              CoordinateTransformation.DMSToDegrees(1, 56, 33.3, false),
                                                              CoordinateTransformation.DMSToDegrees(5, 36, 12.81),
                                                              CoordinateTransformation.DMSToDegrees(1, 12, 7, false));

            bool bType1;
            double separation = AAPlus.AngularSeparation.SmallestCircle(
                CoordinateTransformation.DMSToDegrees(12, 41, 8.63),
                CoordinateTransformation.DMSToDegrees(5, 37, 54.2, false),
                CoordinateTransformation.DMSToDegrees(12, 52, 5.21),
                CoordinateTransformation.DMSToDegrees(4, 22, 26.2, false),
                CoordinateTransformation.DMSToDegrees(12, 39, 28.11),
                CoordinateTransformation.DMSToDegrees(1, 50, 3.7, false), out bType1);

            separation = AAPlus.AngularSeparation.SmallestCircle(CoordinateTransformation.DMSToDegrees(9, 5, 41.44),
                                                                 CoordinateTransformation.DMSToDegrees(18, 30, 30),
                                                                 CoordinateTransformation.DMSToDegrees(9, 9, 29),
                                                                 CoordinateTransformation.DMSToDegrees(17, 43, 56.7),
                                                                 CoordinateTransformation.DMSToDegrees(8, 59, 47.14),
                                                                 CoordinateTransformation.DMSToDegrees(17, 49, 36.8),
                                                                 out bType1);

            Alpha = CoordinateTransformation.DMSToDegrees(2, 44, 11.986);
            double Delta = CoordinateTransformation.DMSToDegrees(49, 13, 42.48);
            Coordinate2D PA = Precession.AdjustPositionUsingUniformProperMotion((2462088.69 - 2451545) / 365.25, Alpha,
                                                                                Delta, 0.03425, -0.0895);

            Coordinate2D Precessed = Precession.PrecessEquatorial(PA.X, PA.Y, 2451545, 2462088.69);

            Alpha = CoordinateTransformation.DMSToDegrees(2, 31, 48.704);
            Delta = CoordinateTransformation.DMSToDegrees(89, 15, 50.72);
            Coordinate2D PA2 = Precession.AdjustPositionUsingUniformProperMotion((2415020.3135 - 2451545) / 365.25,
                                                                                 Alpha, Delta, 0.19877, -0.0152);
            //Coordinate2D Precessed2 = Precession.PrecessEquatorialFK4(PA2.X, PA2.Y, 2451545, 2415020.3135);

            Coordinate2D PM = Precession.EquatorialPMToEcliptic(0, 0, 0, 1, 1, 23);


            Coordinate2D PA3 = Precession.AdjustPositionUsingMotionInSpace(2.64, -7.6, -1000,
                                                                           CoordinateTransformation.DMSToDegrees(6,
                                                                                                                 45,
                                                                                                                 8.871),
                                                                           CoordinateTransformation.DMSToDegrees(
                                                                               16, 42, 57.99, false), -0.03847,
                                                                           -1.2053);
            Coordinate2D PA4 = Precession.AdjustPositionUsingUniformProperMotion(-1000,
                                                                                 CoordinateTransformation.
                                                                                     DMSToDegrees(6, 45, 8.871),
                                                                                 CoordinateTransformation.
                                                                                     DMSToDegrees(16, 42, 57.99,
                                                                                                  false), -0.03847,
                                                                                 -1.2053);

            Coordinate2D PA5 = Precession.AdjustPositionUsingMotionInSpace(2.64, -7.6, -12000,
                                                                           CoordinateTransformation.DMSToDegrees(6,
                                                                                                                 45,
                                                                                                                 8.871),
                                                                           CoordinateTransformation.DMSToDegrees(
                                                                               16, 42, 57.99, false), -0.03847,
                                                                           -1.2053);
            Coordinate2D PA6 = Precession.AdjustPositionUsingUniformProperMotion(-12000,
                                                                                 CoordinateTransformation.
                                                                                     DMSToDegrees(6, 45, 8.871),
                                                                                 CoordinateTransformation.
                                                                                     DMSToDegrees(16, 42, 57.99,
                                                                                                  false), -0.03847,
                                                                                 -1.2053);

            Alpha = CoordinateTransformation.DMSToDegrees(2, 44, 11.986);
            Delta = CoordinateTransformation.DMSToDegrees(49, 13, 42.48);
            Coordinate2D PA7 = Precession.AdjustPositionUsingUniformProperMotion((2462088.69 - 2451545) / 365.25,
                                                                                 Alpha, Delta, 0.03425, -0.0895);
            Coordinate3D EarthVelocity = AAPlus.Aberration.EarthVelocity(2462088.69);
            Coordinate2D Aberration = AAPlus.Aberration.EquatorialAberration(PA7.X, PA7.Y, 2462088.69);
            PA7.X += Aberration.X;
            PA7.Y += Aberration.Y;
            PA7 = Precession.PrecessEquatorial(PA7.X, PA7.Y, 2451545, 2462088.69);

            Obliquity = Nutation.MeanObliquityOfEcliptic(2462088.69);
            NutationInLongitude = Nutation.NutationInLongitude(2462088.69);
            NutationInEcliptic = Nutation.NutationInObliquity(2462088.69);
            double AlphaNutation = Nutation.NutationInRightAscension(PA7.X, PA7.Y, Obliquity, NutationInLongitude,
                                                                     NutationInEcliptic);
            double DeltaNutation = Nutation.NutationInDeclination(PA7.X, PA7.Y, Obliquity, NutationInLongitude,
                                                                  NutationInEcliptic);
            PA7.X += CoordinateTransformation.DMSToDegrees(0, 0, AlphaNutation / 15);
            PA7.Y += CoordinateTransformation.DMSToDegrees(0, 0, DeltaNutation);
        }

        [Test]
        public void RefractionTest()
        {
            double R1 = Refraction.RefractionFromApparent(0.5);
            double R2 = Refraction.RefractionFromTrue(0.5 - R1 + CoordinateTransformation.DMSToDegrees(0, 32, 0));
            double R3 = Refraction.RefractionFromApparent(90);
        }

        [Test]
        public void Parallactic()
        {
            double HourAngle = AAPlus.Parallactic.ParallacticAngle(-3, 10, 20);
            double EclipticLongitude = AAPlus.Parallactic.EclipticLongitudeOnHorizon(5, 23.44, 51);
            double EclipticAngle = AAPlus.Parallactic.AngleBetweenEclipticAndHorizon(5, 23.44, 51);
            double Angle = AAPlus.Parallactic.AngleBetweenNorthCelestialPoleAndNorthPoleOfEcliptic(90, 0, 23.44);
        }

        [Test]
        public void NutationTest()
        {
            Date date = new Date();
            date.Set(1987, 4, 10, 0, 0, 0, true);
            double Obliquity = Nutation.MeanObliquityOfEcliptic(date.Julian());
            double NutationInLongitude = Nutation.NutationInLongitude(date.Julian());
            double NutationInEcliptic = Nutation.NutationInObliquity(date.Julian());
        }

        [Test]
        public void CoordinateTransformationTest()
        {
            Date date = new Date();
            double AST;
            double LongtitudeAsHourAngle;
            double LocalHourAngle;
            Coordinate2D Equatorial;
            Coordinate2D Ecliptic =
                AAPlus.CoordinateTransformation.Equatorial2Ecliptic(AAPlus.CoordinateTransformation.DMSToDegrees(7, 45, 18.946),
                                                             AAPlus.CoordinateTransformation.DMSToDegrees(28, 01, 34.26),
                                                             23.4392911);
            Equatorial = AAPlus.CoordinateTransformation.Ecliptic2Equatorial(Ecliptic.X, Ecliptic.Y, 23.4392911);
            Coordinate2D Galactic =
                AAPlus.CoordinateTransformation.Equatorial2Galactic(AAPlus.CoordinateTransformation.DMSToDegrees(17, 48, 59.74),
                                                             AAPlus.CoordinateTransformation.DMSToDegrees(14, 43, 8.2,
                                                                                                   false));
            Coordinate2D Equatorial2 = AAPlus.CoordinateTransformation.Galactic2Equatorial(Galactic.X, Galactic.Y);
            date.Set(1987, 4, 10, 19, 21, 0, true);
            AST = Sidereal.ApparentGreenwichSiderealTime(date.Julian());
            LongtitudeAsHourAngle =
                AAPlus.CoordinateTransformation.DegreesToHours(AAPlus.CoordinateTransformation.DMSToDegrees(77, 3, 56));
            double Alpha = AAPlus.CoordinateTransformation.DMSToDegrees(23, 9, 16.641);
            LocalHourAngle = AST - LongtitudeAsHourAngle - Alpha;
            Coordinate2D Horizontal = AAPlus.CoordinateTransformation.Equatorial2Horizontal(LocalHourAngle,
                                                                                     AAPlus.CoordinateTransformation.
                                                                                         DMSToDegrees(6, 43, 11.61,
                                                                                                      false),
                                                                                     AAPlus.CoordinateTransformation.
                                                                                         DMSToDegrees(38, 55, 17));
            Coordinate2D Equatorial3 = AAPlus.CoordinateTransformation.Horizontal2Equatorial(Horizontal.X, Horizontal.Y,
                                                                                      AAPlus.CoordinateTransformation.
                                                                                          DMSToDegrees(38, 55, 17));
            double alpha2 = AAPlus.CoordinateTransformation.MapTo0To24Range(AST - Equatorial3.X - LongtitudeAsHourAngle);
        }

        [Test]
        public void SiderialTest()
        {
            Date date = new Date();
            double AST;
            date.Set(1987, 4, 10, 0, 0, 0, true);
            double MST = Sidereal.MeanGreenwichSiderealTime(date.Julian());
            AST = Sidereal.ApparentGreenwichSiderealTime(date.Julian());
            date.Set(1987, 4, 10, 19, 21, 0, true);
            MST = Sidereal.MeanGreenwichSiderealTime(date.Julian());
        }

        [Test]
        public void GlobeTest()
        {
            double rhosintheta = Globe.RhoSinThetaPrime(33.356111, 1706);
            double rhocostheta = Globe.RhoCosThetaPrime(33.356111, 1706);
            double RadiusOfLatitude = Globe.RadiusOfParallelOfLatitude(42);
            double RadiusOfCurvature = Globe.RadiusOfCurvature(42);
            double Distance = Globe.DistanceBetweenPoints(CoordinateTransformation.DMSToDegrees(48, 50, 11),
                                                          CoordinateTransformation.DMSToDegrees(2, 20, 14, false),
                                                          CoordinateTransformation.DMSToDegrees(38, 55, 17),
                                                          CoordinateTransformation.DMSToDegrees(77, 3, 56));


            double Distance1 = Globe.DistanceBetweenPoints(50, 0, 50, 60);
            double Distance2 = Globe.DistanceBetweenPoints(50, 0, 50, 1);
            double Distance3 = Globe.DistanceBetweenPoints(CoordinateTransformation.DMSToDegrees(89, 59, 0), 0,
                                                           CoordinateTransformation.DMSToDegrees(89, 59, 0), 1);
            double Distance4 = Globe.DistanceBetweenPoints(CoordinateTransformation.DMSToDegrees(89, 59, 0), 0,
                                                           CoordinateTransformation.DMSToDegrees(89, 59, 0), 180);
            double Distance5 = Globe.DistanceBetweenPoints(CoordinateTransformation.DMSToDegrees(89, 59, 0), 0,
                                                           CoordinateTransformation.DMSToDegrees(89, 59, 0), 90);
        }

        [Test]
        public void EasterTest()
        {
            Date date = new Date();
            EasterDetails easterDetails = Easter.Calculate(1991, true);
            EasterDetails easterDetails2 = Easter.Calculate(1818, true);
            EasterDetails easterDetails3 = Easter.Calculate(179, false);

            //Test out the AADynamicalTime class
            date.Set(1977, 2, 18, 3, 37, 40, true);
            double DeltaT = DynamicalTime.DeltaT(date.Julian());
            date.Set(333, 2, 6, 6, 0, 0, false);
            DeltaT = DynamicalTime.DeltaT(date.Julian());
        }

        [Test]
        public void DateTest()
        {
            Date date = new Date();
            date.Set(2000, 1, 1, 12, 1, 2.3, true);
            long Year = 0;
            long Month = 0;
            long Day = 0;
            long Hour = 0;
            long Minute = 0;
            double Second = 0;
            date.Get(out Year, out Month, out Day, out Hour, out Minute, out Second);
            long DaysInMonth = date.DaysInMonth();
            long DaysInYear = date.DaysInYear();
            bool bLeap = date.Leap();
            double Julian = date.Julian();
            double FractionalYear = date.FractionalYear();
            double DayOfYear = date.DayOfYear();
            Date.DAY_OF_WEEK dow = date.DayOfWeek();
            Year = date.Year();
            Month = date.Month();
            Day = date.Day();
            Hour = date.Hour();
            Minute = date.Minute();
            Second = date.Second();
            double Julian2 = date.Julian();

            date.Set(1978, 11, 14, 0, 0, 0, true);
            long DayNumber = (long)(date.DayOfYear());
            Date.DayOfYearToDayAndMonth(DayNumber, date.Leap(), out Day, out Month);
            Year = date.Year();
        }

        [Test]
        public void DateSequenceTest()
        {            
            //Do a full round trip test on Date across a nice wide range. Note we should expect
            //some printfs to appear during this test (Specifically a monotonic error for 15 October 1582)

            double prevJulian = -1;

            for (int YYYY = -4712; YYYY < 5000; YYYY++)
            //Change the end YYYY value if you would like to test a longer range
            {
                if ((YYYY % 1000) == 0)
                    Printf("Doing date tests on year %d\n", YYYY);
                for (int MMMM = 1; MMMM <= 12; MMMM++)
                {
                    bool bLeap = Date.IsLeap(YYYY, (YYYY >= 1582));
                    for (int DDDD = 1; DDDD <= Date.DaysInMonth(MMMM, bLeap); DDDD++)
                    {
                        bool bGregorian = Date.AfterPapalReform(YYYY, MMMM, DDDD);
                        Date date = new Date(YYYY, MMMM, DDDD, 12, 0, 0, bGregorian);

                        double currentJulian = date.Julian();
                        
                        Assert.IsFalse((date.Year() != YYYY) || (date.Month() != MMMM) || (date.Day() != DDDD), "Round trip bug with date {0}/{1}/{2} and Julian {3}", YYYY, MMMM, DDDD, currentJulian);

                        Assert.IsFalse((currentJulian != 2299161) && ((prevJulian != -1) && ((prevJulian + 1) != currentJulian)), "Julian Day monotonic bug with date {0}/{1}/{2} and Julian {3}", YYYY, MMMM, DDDD, currentJulian);
                        prevJulian = currentJulian;

                        //Only do round trip tests between the Julian and Gregorian calendars after the papal 
                        //reform. This is because the Date class does not support the propalactic Gregorian 
                        //calendar, while it does fully support the propalactic Julian calendar.
                        if (bGregorian)
                        {
                            CalendarDate GregorianDate = new CalendarDate(Date.JulianToGregorian(YYYY, MMMM, DDDD));
                            CalendarDate JulianDate =
                                new CalendarDate(Date.GregorianToJulian(GregorianDate.Year, GregorianDate.Month,
                                                                GregorianDate.Day));
                            Assert.IsTrue(currentJulian == 2299161 || ((JulianDate.Year == YYYY) && (JulianDate.Month == MMMM) && (JulianDate.Day == DDDD)), "Round trip bug with Julian -> Gregorian Calendar {0}/{1}/{2} with Julian {3}", YYYY, MMMM, DDDD, currentJulian);
                        }
                    }
                }
            }
        }

        [Test]
        public void TopocentricHorizontalPositionMoonTest()
        {
            //Calculate the topocentric horizontal position of the Moon for Palamor Observatory on midnight UTC for the 21st of September 2007
            Coordinate2D Equatorial;
            double Longitude;
            double Latitude;
            double Height;
            double AST;
            double LongtitudeAsHourAngle;
            double LocalHourAngle;
            Date dateMoonCalc = new Date(2007, 9, 21, true);
            double JDMoon = dateMoonCalc.Julian() + DynamicalTime.DeltaT(dateMoonCalc.Julian()) / 86400.0;
            double MoonLong = Moon.EclipticLongitude(JDMoon);
            double MoonLat = Moon.EclipticLatitude(JDMoon);
            Equatorial = CoordinateTransformation.Ecliptic2Equatorial(MoonLong, MoonLat,
                                                                      Nutation.TrueObliquityOfEcliptic(JDMoon));
            double MoonRad = Moon.RadiusVector(JDMoon);
            MoonRad /= 149597870.691; //Convert KM to AU
            Longitude = CoordinateTransformation.DMSToDegrees(116, 51, 45); //West is considered positive
            Latitude = CoordinateTransformation.DMSToDegrees(33, 21, 22);
            Height = 1706;
            Coordinate2D MoonTopo = Parallax.Equatorial2Topocentric(Equatorial.X, Equatorial.Y, MoonRad, Longitude,
                                                                    Latitude, Height, JDMoon);
            AST = Sidereal.ApparentGreenwichSiderealTime(dateMoonCalc.Julian());
            LongtitudeAsHourAngle = CoordinateTransformation.DegreesToHours(Longitude);
            LocalHourAngle = AST - LongtitudeAsHourAngle - MoonTopo.X;
            Coordinate2D MoonHorizontal = CoordinateTransformation.Equatorial2Horizontal(LocalHourAngle, MoonTopo.Y,
                                                                                         Latitude);
            MoonHorizontal.Y += Refraction.RefractionFromTrue(MoonHorizontal.Y, 1013, 10);
            //The result above should be that we have a rising Moon at 17 degrees above the horizon at azimuth 38 degrees east of the southern horizon  
        }

        [Test]
        public void TopocentricHorizontalPositionSunTest()
        {
            //Calculate the topocentric horizontal position of the Sun for Palamor Observatory on midnight UTC for the 21st of September 2007
            Date dateSunCalc = new Date(2007, 9, 21, true);
            double JDSun = dateSunCalc.Julian() + DynamicalTime.DeltaT(dateSunCalc.Julian()) / 86400.0;
            double SunLong = Sun.ApparentEclipticLongitude(JDSun);
            double SunLat = Sun.ApparentEclipticLatitude(JDSun);
            Coordinate2D Equatorial = CoordinateTransformation.Ecliptic2Equatorial(SunLong, SunLat,
                                                                                   Nutation.TrueObliquityOfEcliptic(
                                                                                       JDSun));
            double SunRad = Earth.RadiusVector(JDSun);
            double Longitude = CoordinateTransformation.DMSToDegrees(116, 51, 45); //West is considered positive
            double Latitude = CoordinateTransformation.DMSToDegrees(33, 21, 22);
            double Height = 1706;
            Coordinate2D SunTopo = Parallax.Equatorial2Topocentric(Equatorial.X, Equatorial.Y, SunRad, Longitude,
                                                                   Latitude, Height, JDSun);
            double AST = Sidereal.ApparentGreenwichSiderealTime(dateSunCalc.Julian());
            double LongtitudeAsHourAngle = CoordinateTransformation.DegreesToHours(Longitude);
            double LocalHourAngle = AST - LongtitudeAsHourAngle - SunTopo.X;
            Coordinate2D SunHorizontal = CoordinateTransformation.Equatorial2Horizontal(LocalHourAngle, SunTopo.Y,
                                                                                        Latitude);
            SunHorizontal.Y += Refraction.RefractionFromTrue(SunHorizontal.Y, 1013, 10);

            //The result above should be that we have a setting Sun at 21 degrees above the horizon at azimuth 14 degrees south of the westerly horizon  
        }

        private void Printf(string message, params object[] args)
        {
            Utility.printf(message, args);
        }
    }
}
